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Abstract 



We find the static displacement, stress, strain and the modified Coulomb failure stress pro- 
duced in an elastic medium by a finite size rectangular fault after its dislocation with uniform 
stress drop but a non uniform dislocation on the source. The time-dependent rate of triggered 
earthquakes is estimated by a rate-state model applied to a uniformly distributed population of 
faults whose equilibrium is perturbated by a stress change caused only by the first dislocation. 
The rate of triggered events in our simulations is exponentially proportional to the Coulomb 
stress change, but the time at which the maximum rate begins to decrease is variable from frac- 
tions of hour for positive stress changes of the order of some MPa, up to more than a year for 
smaller stress changes. As a consequence, the final number of triggered events is proportional to 
the Coulomb stress change. The model predicts that the total number of events triggered on a 
plane containing the fault is proportional to the 2/3 power of the seismic moment. Indeed, the 
total number of aftershocks produced on the fault plane scales in magnitude as 10 M . Including 
the negative contribution of the stress drop inside the source, we observe that the number of 
events inhibited on the fault is, at long term, nearly identical to the number of those induced 
outside, representing a sort of conservative natural rule. Considering its behaviour in time, 
our model doesn't completely match the popular Omori law; in fact it has been shown that 
the seismicity induced closely to the fault edges is intense but of short duration, while that 
expected at large distances (up to some tens times the fault dimensions) exhibits a much slower 
decay. 
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1 Introduction 

Studying the space-time interaction of earthquakes is important not only for the comprehension 
of the phenomenon of earthquakes, but also for its possible application in the mitigation of 
earthquake risk. 

In this paper we aim to the physical modelling of the interaction between seismic events, 
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seeking a relationship between the source parameters of an earthquake and the rate of all other 
earthquakes that follow it, a phenomenon empirically and variably labelled with the names of 
seismic acceleration, foreshock and aftershock sequences, swarms, earthquake clustering, etc. 
A purely statistical approach to this p henomenon exists with the name of "epidemic model" 



(Ogata. 



Console et al 



20031 ). In the epidemic model each earthquake source is 



supposed capable of increasing the probability of new earthquakes according to a kernel space 
and time probability distribution decreasing from the source. The epidemic model is applicable 
not only to the aftershock phenomenon but also to the foreshock's. However, despite it provides 
an accurate statistical rapresentation of earthquake interaction, the physical interpretation is 
ignored. 

In order to formulate a more physical description of the clustered seismicity, we consider 
that seismic events modify the stress field around the causative fault. Our main goal in this 
work is to observe the spatio-temporal variation of the seismici ty rate after an event intro- 



ducing in t he algorithm a ra te-state model approach derived by 



(toieterich 



1992, 



Ruinal (|l983h and Dieterich 



1994). In particular, we seek a physically based relationship between 



the magnitude of the causative fault and the number of its aftershocks. We aim to create 
a link be tween a physi c al mo delling of the stress transfer and the epidemic model approach 



(see also 



Console et al 



( 20041 ) ). Some recent studies put in evidence that sudden stress vari- 
ations, even of small magnitude, can produce large variations of the seismicity rate. It is 
recognized as a phenomenon of triggering, rather than induction, on faults that are already 



subject to relevant stress. The seismic rate increases in general where t 
Coulomb Stress Change) is positiye, according to the Coulomb mode 



le stress change (called 



1988; 



Boatwright and CoccoL 



1996 



Stein et al 



1997 



Gomberg et al. 







Mcndoza and Hartzell, 

n — T 

0; Harris, 



1998, 



200( 



1998 



Kin g and Cocco 



2001 



Stein. 



1999; 



Tod a and Stein 



200 



Kilb et al 



2002; 



Belardinelli et al.1 



2003r ?). These studies were able to give a physical interpretation for earthquake interaction 
observed in specific real cases, but were not suitable for the general application of the model in 
predictive way as it is necessary for the information contained in a seismic catalog. In partic- 
ular, the Coulomb stress change criterion doesn't allow to model the hyperbolic time decay of 
the aftershock rate after a mainshock, known as the Omori Law. In this simple model all the 
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next earthquakes would be advanced by a time wich depends on the stress change and on the 
loading rate but doesn't depend on the physical properties of the fault. The experience shows 
that the Earth doesn't rather behave in this way. The most popular empirical description of 
the phenomenon is the Omori law, stating that the the aftershock rate decays in time as f 1 . 
Among the various theories modelling the Omori law, we have taken in consideration only the 
one that involves a particular constitutive behaviour of faults. Therefore, we have not consid- 
ered, even without rejecting them in principle, other hypotheses that predict a time variation 
in stress, like the viscoelasticity and the diffusion of fluids in the medium. The rate-state model 
for earthquake nucleation seems capable of substantially explaining all the phenomenology. In 
this model slip doesn't occur instantaneously at the exceeding of a threshold stress as for the 
Coulomb-Amonton criterion, but follows a more complex time history with different phases. 
This time history depends on the physical properties of the fault too. ? have shown that the 
rate-state Dieterich model may be derived from simple general i deas. The applica tion of the 
rate-and-state model to the evolution of seismicity is not new, e.g., 



To da et al 



( 19981 ) have used 



this model to predict the spatio-temporal seismicity rate after the Kobe earth quake. Other ex- 



ampl e s of applica tion of the same model can also be found in the literature ([Toda and Stein . 



200(1 



2002 



2003). Following the rate-state model it has recently been possible to simulate 



seismicity quite realistically, accounting for the rate of eyents induced by stress changes even 



2003 



Ziv and Rubin . 



2003). 



at large distance from the inducing earthquake (JZiyj, 

Basing our work on the classical theory of the elasticity, we start from the formulation of 
the stre ss released by a rectangular fault with uniform stress-drop; then we apply the rate-state 



model ( Dieterich 



19941 ) obtaining the complete time and space distribution of the seismicity 



induced by a given fault. It is important to remark that for present purposes we neglect the 
interaction between subsequent events, supposing that only the first event perturbs the stress 
field. We also neglect the effect of the free surface of the Earth. In this study we limit our 
interest to the behaviour of induced seismicity in space and time and its relationship with the 
fault's parameters. Because of the qualitative natu re of thi s work we don't use the expressions 



derived, in a similar but more complex context, by 



Okadal (jl985a[ ). The results of simulations 



are analysed to find out the scaling relationships existing between the free parameters of the 
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model and the expected seismicity in a way that will allow the validation of the model by real 
observations. 



2 Elastic model for a rectangular fault 



The theoretical approach of this work has been made as simple as possible, still preserving the 
capability of modelling the stress transfer from a seismic source to another. 

Suppose that a fault is embedded in a homogeneous, isotropic, elastic medium where the 
stress tensor a is uniform. At a particular moment, the fault slips generating an earthquake, 
and suppose that the earthquake fully releases the component of the traction parallel to the slip 
vector across the fault. We call Act (the stress drop) the uniform negative change of traction 
on the fault. It is well known that the slip distribution on a fault in similar conditions is not 
uniform. The slip distribution will satisfy the theory of the elasticity at the new equilibrium, 
and it will be zero on the edges of the fault. It is also well known that the shear stress on the 
fault plane outside the edges of the fault increases significantly. 

The analytical solution for the slip distribution on a fault at the equilibrium doesn't exist 



except that for very sim ple geometries, such as a rectangular f ault o 



by 



le ngth (|K nop off , 



1959; 



constant width and infinite 



Udas 



1999). As mentioned 



1957) and for a circular fault ([Keilis-Borok 
Kostrov and DaJ ( 1998 ). the analytical solution is not known even for a geometry as simple 



as a rectangular shape of width W and length L. In the literature it is possible to find the 
computation of the stress field generated by an ideal rectangular fault with uniform slip in 



an infi nite homogeneous space or semispace (jChinnery 
1985hl l 



1961 



Twasaki and Sato 



12Z9J 



Okada. 



However, as mentioned above, this is not a situation of equilibrium fulfilling all the 
physical constraints imposed by the theory of the elasticity. For this reason, here we make 
use of a non uniform slip distribution that is approximatively compatible with a uniform stress 
drop on the fault. 

To calculate the static displacement, stress and strain, we start from the components % of 
the static displacement vector u\ (x) produced by an elementary force directed as the j axis in 
any point of an infinite, elastic, homogeneous and isotropic medium, known as the Somigliana 
Tensor. From the expressions valid for a single force source, those for a couple of forces of given 



moment are derived by differentiation with respect to the direction of the arm of the couple of 
forces. Then the expressions for a source represented by a double couple of forces are obtained 
by summing a pair of two co-planar couples orthogonal to each other. Once the expressions 
for the displacement are known, the stress cr^-, the strain and the tractions are computed 
by their definitions. Indeed we also can compute the modified Coulomb failure stress from the 
relation: 

a c = T s -/j'T n , (1) 

where T s and T n are respectively the traction parallel to the slip vector and the traction normal 
to the fault plane (negative for compression); // is the effective coefficient of friction incorpo- 
rating the pore fluid pression that is commonly assumed to be proportional to normal stress 
change. The effective coefficient of friction in coseismic stress change studies has varied from 
0.0 to 0.75, with an average value of y! = 0.4 often used. This is the value we have used in our 
simulations. 

Our model of an extended source is drawn euristically from that of a rectangular fault 
infinitely extending in one direction and that of a circular fault, and it is characterized by a 
slip distribution as: 



Au = — 



\ 



(f) : 



(I) 2 - v 2 



LW 
4 



(2) 



defined for — W/2 < x < W/2 e — L/2 < y < L/2, where x, y are the coordinates on the fault 
plane whose origin is coincident with the center of the rectangle, Act is the stress drop and fi the 
shear modulus (rigidity) of the elastic medium. Equation (j2J) satisfies the condition of zero slip 
on the edges of the rectangle. The slip distribution is of a form very similar to that obtained 
numer ically by discretization, except for the cases when one dimension is much larger than the 



other ( Kostrov and Da; 



1998|) • Its analytical form is similar to the solution for a circular fault. 
Moreover, for a fault of square shape (W = L) it achieves a seismic moment consistent with 
those obtai ned analitically for two ci rcular faults, respectively inscribed and circumscribed to 



the square. 



Bonafede and Neri 



(J2000J) have shown that, when imposing a unidirectional traction 
release in the strike or dip direction, a minor component of slip is present, over the fault plane, 
even in the direction perpendicular to the released traction. In consideration of the qualitative 



character of this study, we have neglected this component. 



For our numerical applications we discretize the continuous slip distribution of equation 
P| with a set of point sources densely distributed on the fault. The more numerous are the 
point sources, the more accurate will be the approximation. Each element (p, q) of the fault is 
characterized by its scalar seismic moment m pq that, by the theorem of the representation is 
given by: 

m pq = H-^ U pq^pqi (3) 

where Au pq is the slip on the element (p, q) and S p5 is its elementary area. 

If u pg (x) is the displacement produced by the slip on the element (p, q) in the point x, the 
displacement resulting from the combination of all the elements is computed by: 

M N 

uW = £5]u M (x), (4) 

p=l q=l 

where u pg (x) is computed by the relation valid for a point source of moment m pq . 
Following equation (j2J) the slip Au pq is written as: 



Au nn = Au r 



(2) X lp (T) X lq 



~4~ 



(5) 



where, in this and xz q are the x\ and x 3 coordinates of the element (p, q) respectively, 

and: 

Au rn ax = ~ ^ ^ (6) 

is computed through the other three fault parameters. These parameters are also connected to 
the scalar moment M through the relations 

2 2 

M = ^WL = \-uAu max WL = l~Aa{WLf' 2 (7) 
lo 62 

where Au is the average slip on the fault. 

The above mentioned theoretical approach has been implemented in a computer code which 
allows the user to chose among several options: - The focal mechanism characterizing the 
causative fault and the induced events constrained by the strike 0, the dip 5 and the rake A 
angles; - Dimensions and stress drop, seismic moment, or magnitude of the fault; - Type of 
output: Displacement (3 components), stress or strain (6 independent components), traction 



parallel to the slip vector and perpendicular to the fault plane (2 component) and Coulomb 
failure stress (1 component). 

The output values are given on a rectangular grid of points, whose orientation (horizontal, 
or vertical parallel to one of the co-ordinate axes), dimensions, and spacing are also selected by 
the user. See Figure 1 for the conventions we adopted for the reference system and the angles. 

In this work we make, for sake of simplicity, the assumption that all the aftershocks have 
the same focal mechanism of the causative fault. 

Figures 2 and 3 show the component of the stress field parallel to the slip <j\i and the 
normal traction T n produced by a vertical strike slip fault whose parameters are listed in Table 
H Note that these figures show the stress change and not its absolute value. 

Figure 2 clearly shows that the shear stress a\2 is maximum near the edges of the fault, in 
particular along the slip direction x\, and decreases with increasing the distance from the fault. 
Inside the fault the stress change is negative and approximately equal to Act. 

Figures 4, 5 and 6 show respectively au, T n and the Coulomb failure stress a c produced by 
a vertical strike slip fault on the horizontal plane. 



3 The rate-state model 



As already stated, we assume that the time-dependent behaviour of the seismicity triggered by a 



Co ulomb stress change in a pop ulation of 



by 



Ruinal ()1983f ) and Dieterich ((Dieterich 



aults is described by the rate-st ate model i ntrod uced 



1986 



1992, 



19941 1 . According to 



Dieterich (1994 the 



rate R of earthquakes due to a At Coulomb failure stress is given by: 

R = - — T^O, 



(8) 



whose various parameters are defined in Table El Note that equation (JBJ) doesn't predict the 
magnitude of the earthquakes, for the distribution of which we assume the classical Gutenberg- 
Richter frequency-magnitude relationship. However, both R and r refer to the number of 
earthquakes in a specific magnitude range. 
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The modelling of the space-time distribution of seismicity triggered by an earthquake is 
straightforward, introducing the expression for the Coulomb failure stress a c in equation (jHJ 
for At. We make the assumption f = f T . The meaning of this assumption is that the primary 
earthquake doesn't change the long-term stress rate due to the tectonic driving forces in the 
environment. 

It is easily understandable from (jHJ) that the integral of R over infinite time diverges. 
This is due to the fact that the limit of R for t — > oo is the background rate r. We prefer 
to deal with the net triggered seismicity rate R' = R — r. The plot of R' versus time for 



Dieterichl (J1994)). It shows that the initial 



different values of A r is shown in Figure 7 (see also 
value of the triggered seismicity rate, immediately after the triggering event, is proportional to 
the exponential of the stress variation but the time period over which R(t) is maximum and 
relatively constant decreases exponentially by the same fa ctor for increasi ng Coulomb stress 



Dieterichl (J199J), can be obtained 



change. This time period, in substantial agreement with 
computing the intersection between R'(t = 0) and the stra ight line corresponding to the limit 



of R' when Ai 



oo. The difference with 



Dieterich 




as Figure 8 shows, is that we have 



preferred to make use of R', rather than R, and we ignore the case f = which we judge not 
realistic. Such time period ranges from few hours for a variation At = 5MPa to nearly one 
year for At = IMPa. Figure 7 finally shows that the total duration of the triggered seismicity 
is the same but the shape of the decay is different for different At. 

The integration of R' over infinite time leads to the not completely trivial result thet the 
total number of tr i ggere d events in an area of constant stress change is proportional to At 



( Helmstetter et al 



2004): 



R'(t)dt = ^At = -At. 
w Aa t 



(9) 



by 



This relationship predicted by the friction model had been observed for Landers Aftershocks 



Gross and Kisslingerl (|1997l) 



Consider the integral over time of the triggered seismicity rate R', once for a positive At, 
and then for a negative At of the same size, and we sum the two functions. We obtain a 
quantity that intially has a positive increase and then tends back to zero (Figure 9). This 
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result can be interpreted as a sort of conservative law of the seismicity, by saying that the total 
number of aftershocks generated in a place close to a fault can be balanced, at long range, by 
the total number of those events that are inhibited by the stress drop in a place internal to the 
fault if the stress drops are equal in absolute value. It seems to us that this circumstance has 
been ignored, so far, in all previous studies on this subject. 

For numerical applications it is necessary to define the value of the various parameters 
appearing in equation (JSJ). Table El contains a list of values inferred from several geophysical 
considerations and used in our simulations. The reference seismicity rate, r, is a quantity that 
should be known experimentally for a given area, and it is related to the average strain rate f of 
that area. We use for r a uniform value of four events of magnitude exceeding 3.0 per year per 
1000 km 2 , and for f a value of 5 KPa per year. These are simplifications of the reality where r 
and f are not geographically uniform. These parameters seem however reasonable for an area 
of moderate seismicity. The reason why r is expressed as a number of events per unit area and 



time is given below. The parameter A of the co nstitutive 



aw h as a value ranging from 0.005 



19951 ). Generally, simulations that 



to 0.015 obtained from laboratory experiments (jDieterich . 
most nearly resemble earthquakes in nature were obtained with rather small values of A, near 
to A fa 0.001. The value of Aa has been evaluated for different sequences of earthquakes by 



several autors (jBelardinelli et al 



2003 



Stein 



1999). The parameter a is the normal pressure 



on the fault (of the order of some tens of MPa). The characteristic time of aftershocks, t a , 
depending on the other parameters through the relation t a = (Aa)/f, in this context is of the 
order of several tens of years. This value is essential to model a long term triggered seismicity, 
because the effect of the stress change disappears completely after a time of the order of the 
double of t a . 



4 Application to a rectangular fault 



Dieterichl ([19941 ) discussed in detail the case of a finite circular dislocation with uniform stress 
drop. In our simulations we have taken into consideration a vertical rectangular fault with 
horizontal slip in the x\ direction, embedded in a homogeneous, isotropic, elastic medium. The 
values assumed for the fault parameters are again those reported in Table H Assuming a con- 
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stant stress drop Act, equation (0) shows that the scalar seismic moment Mo is proportional to 
the quantity (WL) 3 ^ 2 , that can be considered a sort of equivalent source volume. Scaling W and 
L by a given factor, M scales accordingly, and also the shape of the stress distribution in space 
scales by the same factor. Assuming that the background seismicity rate is given as the number 
of events per unit volume and time, and integrating over the space, we find that the number of 
events triggered by the earthquake is proportional to the cube of the linear dimensions, i.e. to 
the seismic moment M , or to 10^ M . In this study we make the different hypothesis that the 
aftershocks of an earthquake occur mainly on the fault plane. This hypothesis comes from the 
consideration that the fault plane is subject to the maximum sheare stress and that it is a sort of 
plane of weakness where the triggered seismicity is more likely to occur. In this case, the num- 
ber of triggered events would be proportional to M^ 3 , or to 10 M . If it was true, and neglecting 
other kinds of factors, one earthquake of given magnitude would have, all togheteher, a number 
of aftershocks ten times larger than one earthquake of one magnitude unit smaller. Assuming, 
then, that the b value of the Gutenberg- Richter frequency- magnitude relation is equal to 1, it 
would follow that the earthquakes of a magnitude class produce the same number of aftershocks 
as the earthquakes of the magnitude class one unit smaller. This hypothesis is supported by 



the statistic al analysis of real observations (IFelzer et al 



Htsii 



2001 



Yamanaka and Shimazakil . 



1990; 



1969) and will be used in the following of the paper. 
A first group of simulations concerning the seismicity triggered by a fault as a whole have 
been performed to study the variation of the total number of triggered events (in space and 
time), Ntot, versus the main free source parameters: the stress drop Act, the linear dimensions L 
e W and the magnitude M. We assume that the magnitude M is proportional to the logarithm 



Kanamori and Anderson 



(1975) 



of the seismic moment M (measured in Nm) according to the 
relation: 

log M = 9.1 + 1.5M. (10) 

Tables 01 El and El show the results obtained assigning a fixed value for one of the parameters 
at a time, and letting the others change. These simulations have been carried out by summation 
of the single contributions of the elementary cells constituting the area surrounding the fault. 
Only the seismicity triggered on the fault plane, but outside the fault edges, has been considered. 
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Looking at Table0]one can notice that, under our hypoteses, for a constant stress drop, the total 
number of triggered events increases about proportionally to the fault area LW. Theoretically 
it should be exactly proportional to the area of the source, but in the calculations there are 
numerical imprecisions due to the grid size approximations. These numerical simulations show, 
also, that the number of aftershocks is proportional to the 2/3 power of the scalar seismic 
moment or to the exponential of the magnitude (10 M ), in agreement with the above theoretical 
arguments. Table 03 clearly shows the proportionality existing between the stress drop Aa and 
N to ti having fixed the linear dimensions of the fault. Furthermore, the results reported in Table 
Elgive credit to the hypothesis that, for a constant magnitude, the earthquakes of smaller linear 
dimensions generate a larger number of triggered events. 

Let's mention, lastly, the results of a simple numerical test. We just have theoretically 
shown in the previous section that the total number of events triggered by a stress change 
At per unit area is, after a time much larger than the characteristic time t a , identical to the 
number of events inhibited by a negative Coulomb stress change — At on the same area. The 
numerical simulations about an extended source reflect this idea in a more general way. In fact, 
if we take the total number of aftershocks produced in a time interval of 100 years by a fault 
of 100 km 2 on the portion of an area of 10 4 km 2 external to it, we find 353.2 events, while the 
same computation made including the negative contribution of the area internal to fault itself 
gives only 15.6 (less than 5%). The idea of an overall null balance is, then, substantially met: 
the phenomenon of clustering doesn't change the long term seismicity rate of a large area. This 
is in contrast with the opinion that trigg er zones produce many mo re shock s than are missing 



from the shadows ( Toda and Stein 



2003). We neverthless regard for 



Marsanl (J2003J) that put in 



evidence the practical difficulty to observe in nature the lack of events (the quiescence) mainly 
when looking at weakly active regions and short timescales. 



5 Temporal behaviour of the triggered seismicity 



Let's now consider in more detail t he dipendence on time of the n umber of 



predicted by the present hypothesis. 



Kagan and Jackson! (jl99lh and 



trigge red events 



Dieterichl (jl994f ) found, in 



their empirical works, that shallow aftershocks, above a given magnitude cutoff, decay globally 
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(for different distance intervals) in about 10 years and this time decay appears to sistematically 
decrease with increasing depth. In our model the behaviour in different zones depends on the 
different values of the stress changes, that decrease with the increasing distance from the edges 
of the fault. It can be then noted that the plateau of the plot of the seismicity rate in time 
(see Figure 7) has a duration that increases with the increasing distance from the source. This 
result depends on our assumption of considering only aftershocks on the mainshock fault plane. 
It can also be biased by the fact that we consider only the effect of static stress step At and we 
ignore the interaction between events. In spite of that, we believe that our model can decribe 
a behavior not so far from reality. Therefore we keep as our reference the model of fault given 
in Tabled again. 

As shown in Figure 10a, in the close neighbouring of the fault (within 0.5 km from the 
perimeter of the fault), most of the aftershocks occur within one day after the occurrence of 
the inducing earthquake, and more than 1/3 of them occur already in the first hour. This is 
clearly not the case for further distance from the source: in the slice between 0.5 and 8.5 km 
from the fault, about half of the activity is exploited within about one year. The further we 
go apart from the source, the longer is the time interval necessary for the exploitment of the 
triggered activity. In fact, in the distant zone (from 20 to 100 km from the fault) most of 
the triggered events occur longer than 20 years after the inducing earthquake (Figure 10b) . 
The behaviour shown in Figure 10a and 10b is illustrated also by the snapshots of the spatial 
distribution of the seismicity rate at different times. Figures 12, 13 and 14 clearly show how 
the distribution of aftershocks resembles that of the Coulomb stress change (Figure 2) and how 
the time decay is different in the different slices. These also show how shadow zones (where the 
seismicity is inhibited) have an effect longer in time. It should be noted again that we make 
the assumption that the aftershocks occur mainly on the fault plane and that these have the 
same focal mechanism of the causative event. 

We are now interested to check if the generalized Omori law is suitable to describe the 
temporal behaviour of R'(t) predicted by our model in every spatial slice and, in such case, 
by means of which free parameters. We refer to the gene ralized Omori law as to the formula 



describing the decay of aftershock rate after a mainshock (|Utsu et al 



19951 ) through three free 
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parameters a, c e p: 



R'(t) 



(c+ty 

For this purpose we make use of an algorithm for the least squares best-fit of the sets of data 
reported in Figure 10 by the integral of (jllj) over the time: 



R\t')dt' = a 



1 — p 
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The results obtained for the best-fit are shown in Tabled It is evident that the values of the 
parameters obtained are quite different from case to case (for instance, the c parameter is ex- 
tremely small for the slices closest to the fault edges), and that the standard deviation is larger 
for the slices that include the wider range of distances from the fault. These results show that 
the same set of parameters in equation (|12jl can not fit the synthetic values obtained for different 



zones. The observation that the Omori p para meter varies spatially in rea 
has been documented in a number of papers (IWiemer an d Katsu matal. 



2002; 



Wiemer 



2000; 



Wiemer et al 



2002; 



Wiemer and Wyss 



2002). 



aftershocks sequences 



1999; 



Enescu and Ito 



Wiemer et al 



I 



(2002$ re- 



ported a case of significantly d ifferent c values in th e northern and southern Hector Mine 



aftershock volume, respectively. 
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Enescu and ltd (J2002J) in their study of aftershock activity of 



the 2000 Western Tottory earthquake found c- values very close to except for larger magnitude 
earthquakes, arguing that this feature could result from incompleteness of data, or might also 
reflect the complexity of the rupture process. 

The situation is illustrated also by the plots of Figures 15 and 16, showing the graphical 
comparison between the synthetic data and the relative theoretical approximation for different 
slices. 



6 Conclusions 

In the development of this model some very crude approximations have been made. Two of the 
most important are (a) having neglected the dishomogeneneity of the stress drop on the fault, 
and consequently the very complicated pattern of slip (this must certainly have influence on 
the aftershock pattern: some aftershocks may even occur inside the fault, where the gradient 
of slip is high); (b) having neglected the interaction between subsequent events. 

14 



The first step of our modelling, concerning the stress change in the medium fol 
earthquake, has given re s ults c o mparable to those of prev ious papers as, for instance, 



(1961). 



King and Coccol 1)200 lh . 



Belardinelli et al. 



owing an 



Chinnery 



( 20031 ). This substantial agreement, though 



the simplifications introduced in this specific case of rectangular source, supports the validity 
of our methodology. 

As to the i nteraction among seismic events, the application of the rate-state model of 



Dieterichl (J199J) has shown how the total number of triggered events per unit area is propor- 



tional to the Coulomb stress change. Nevertheless, the time constant by which the rate of 
events decreases has a large variability, ranging from a fraction of hour for a Coulomb stress 
change of several MPa, up to tens of years for stress variations smaller than 1 MPa. 

A first set of simulations carried out with a model of rectangular fault, based on the 
hypothesis that most of the aftershocks occur on the same plane of the fault, has shed light 
on the relationship between the main source parameters and the number of triggered events. 
The simulations lead to the conclusion that many earthquakes of small magnitude can produce, 
in total, a number of aftershocks comparable to that of fewer larger earthquakes. Our model 
predicts, in fact, that, for a constant stress drop, the number of triggered events is roughly 
proportional to the seismic moment at the power of 2/3, or equivalently to 10 M . 

An important aspect of the model concerns the time-space behaviour of the triggered seis- 
micity. The s imulations show that the application of the rate-state model and equation (jH)) 



introduced by 



Dieterichl (|l994f ) to the static model of stress change has as a consequence some- 



thing more complicate than the conventional Omori law for the temporal decay of aftershock 
rate. In fact, the seismicity is expected to be intense (in accordance to the stress change,) and 
have a maximum constant value, close to the fault edges of the primary earthquake, but the 
time necessary for the starting of decay is shorter. On the contrary, the seismicity expected 
at larger distances (up to 10 times the linear dimensions of the fault) is weak, but the time 
necessary for the decay from the maximum value is longer (even some tens of years). We can 
conclude that the temporal behaviour of aftershock rate varies in space and this is theoretically 
because " c" increases with the distance from the fault edge. Moreover, the model predicts that 
the total long term production of aftershocks, including wider areas, is not negligible respect 
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to the short term production. In this respect, the regional background seismicity could be 
interpreted as a sort of noise, or memory, due to the superposition of the effect of many older 
earthquakes. The time decay described by the popular Omori law could be interpreted as an 
apparent average result of the contribution from the various areas on the plane containing the 
primary fault. 

Let's list, in short, the most important results concerning the static effect of a dislocation, 
as following: 

• The number of triggered events is proportional to the area of the fault, for a constant 
stress drop; 

• the number of aftershocks is proportional to the 2/3 power of the static seismic moment, 
for a constant stress drop; 

• the logarithm of the total number of aftershocks scales linearly with the magnitude, for 
a constant stress drop; 

• the number of triggered events is proportional to the stress drop on the fault, if its linear 
dimensions are kept constant. 



The rate-state model ( Dieterich 



19941 ) used in our algorithm for the time distribution of the 
seismicity produced by a dislocation on a rectangular fault, has allowed a partial justification 
of the Omori law. We outline the following main conclusions obtained from the simulations: 

• The seismic activity is very intense during the first few hours or days after the occurrence 
time of the primary earthquake, in the rectangular slice close to the edges of the fault; 

• The immediate decay of this intense activity is followed by a period of time during which 
the activity at larger distance (of the order of the linear dimensions of the fault) is nearly 
constant. For longer time (several years) the total number of triggered events in the 
external zone becomes comparable with that of the most internal one. 

• the number of events triggered over a long time scale at distances larger than the fault 
linear dimensions is significantly large; 
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• the seismic activity inside the fault drops to negligible values at the occurrence of the 
primary earthquake and returns to normal values only after a time much longer than the 
characteristic time (several tens of years in our simulations); for a long time range the 
number of events inhibited inside the fault is comparable to that of the events triggered 
outside: the net balance can be considered null; the seismicity rate changes in space and 
time, but it doesn't have influence on the rate averaged over very large intervals of space 
and time. 

In light of the consequences derived from our model, both positive and negative perturba- 
tions of the seismic activities resulting from the interactions among earthquakes would last for 
a time of the order of decades. In this respect, it seems unrealistic to define and observe what 
is commonly called " 'background seismicity" ' of an area. 

We can conclude that the model here analysed, though its very simple conception, gives a 
physical justification both to the very popular phenomenon of short term earthquake clustering 
(aftershocks and foreshocks in strict sense) and to that of the long term induction, also observed 
in various occasions. 

We are confident that the refinement of the model and its validations with real catalogs 
will bring to even more interesting consequences. 
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Tables 



Table 1: Source parameters used in numerical applications. 



Parameter 


Value 


Strike, dip and rake 


0,90,0 

7 7 


Dimension of the fault 


10 kmx6 km 


Spacing on the fault 


0.1 km 


Dimension of the grid 


30 A;mx30 km 


Spacing on the grid 


0.5 km 


Magnitude 


6 



Table 2: Constitutive parameters used in Dieterich's relation. 



Symbol 


Description 


r 


reference seismicity rate 


f 


shear stressing rate 


f r 


reference stressing rate 


At 


earthquake Coulomb stress change 


A 


fault constitutive parameter 


a 


normal stress 


t 


time 


ta 


characteristic time for seismicity to return to the steady state 
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Table 3: Values of the parameters used in the numerical applications of the model for induced 
seismicity. 



Parameter 


Value 


r 


4 events /(y ■ 1000/cm 2 ) 


f 


5 KPa 




5 KPa 


A 


0.008 


a 


30 MPa 


t a 


48 y 
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Table 4: Relation between fault parameters and aftershock production. The value of Act is 
fixed at 4 MPa; the grid used in this computation has an area of 10 4 km 2 , the spacing on it is 
0.5 km. The spacing on the fault is 0.05 km. 



L ■ W [km 2 ] 


M 


M 


[Nm] 


Ntot 


4 


4.6 


9.9 


•10 15 


13.4 


25 


5.4 


1.5 


•10 17 


78.2 


100 


6.0 


1.2 


•10 18 


263.2 


225 


6.4 


4.2 


•10 18 


532.8 



Table 5: Relation between fault parameters and aftershock production. The area of the fault 
is fixed at the value of 100 km 2 , the spacing is of 0.05 km; the grid used in the computation 
has an area of 17 x 17 km 2 , the spacing on it is of 0.5 km. 



Act [MPa] 


M 


M 


[Nm] 


Ntot 


1 


5.6 


3.1 


■ 10 17 


51.8 


2 


5.4 


6.2 


•10 17 


103.4 


3 


5.9 


9.3 


•10 17 


155.2 


4 


6.0 


1.2 


•10 18 


206.8 


5 


6.1 


1.6 


•10 18 


258.6 



Table 6: Relation between fault parameters and aftershock production. The value of M is fixed 
at 6.0; the output grid is of 20 kmx20 km, the spacing is 0.5 km. The spacing on the fault is 
0.05 km. 



L ■ W [km 2 ] 


Act [MPa] 


Not 


4 


510.7 


overflow 


25 


32.7 


597.6 


100 


4.1 


223.2 


225 


1.2 


107.8 
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Table 7: Results obtained for the best-fit of the sets of data reported in Figure 10 by the 
relation (jl'2j) . The symbol rms indicates the standard deviation. 



Slice 


a [y?- 1 } 




V 


rms 


Co- 


'- 0.5) km 


3.24 ±0.17 


0.37- 10~ 7 ± 0.2- 10~ 6 


1.02 ±0.02 


1.86 


co - 


'- 8. 5) km 


11.0 ±0.76 


0.23- 10~ 7 ± 0.2- 10~ 6 


0.91 ±0.032 


9.16 


(0.5 


-r- 8.5)km 


83.7 ±70.0 


3.1 ± 1.96 


1.44 ±0.19 


3.3 


(o- 


± 20)km 


34.3 ±8.0 


0.13- 10-^0.019 


1.14 ±0.077 


9.22 


(o- 


- 100)km 


36.8 ±9.1 


0.15- 10-^0.023 


1.12 ±0.08 


10.8 


(20- 


± 100)km 


12.5 • 10 3 ± 10.9 • 10 3 


61.0 ±5.43 


2.8 ±0.96- 10~ 3 


0.45 
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Figure captions 

Figure (1): Conventions adopted for the reference system and for the parameters of the 
orientation of a fault: the strike 0, dip 5 and rake A angles and the unit vectors n and 1. The 
x 3 axis is positive towards the Nadir direction. 

Figure (2): Shear stress, <7i 2 , parallel to the slip direction for a strike-slip fault projected 
on the fault plane (xix$, x 2 = 0.2 km). 

Figure (3): Tractions normal to the fault plane, T n , projected on the plane X1X3 (the same 
plane of the fault; (x 2 = 0.2 km). 

Figure (4): Shear stress, a± 2 , parallel to the slip direction for a strike-slip fault on a plane 
perpendicular to the fault (0:1X2). 

Figure (5): Tractions normal to the fault plane, T n , projected on the plane perpendicular 
to the fault (2:1X2) . 

Figure (6): Coulomb stress change, At, projected on the plane perpendicular to the fault 
(xix 2 ). It is a combination of <7i 2 and T n on the same plane. 

Figure (7): Distribution of the triggered seismicity against time for different values of 
induced shear stress. The values of the parameters used are: r = 2 events /y, At = 3 MPa, 
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Ao = 0.24 and f = 5 KPa/y. 



Figure (8): Comparis on between the behaviour of the rate against time in the two cases 



treated by 



Dieterichl (jl994l ) and the case that we have selected in our study. 



Figure (9): Distribution of the number of expected events in four different cases: (a) for 
a positive value of earthquake shear stress equal to 5 MPa, (b) for a negative value of it equal 
to -5 MPa, (c) the sum of (a) and (b) and (d) the integral of (c) over time. The values of r and 
t a are fixed at 1. 



Figure (10): Rappresentation of the total number of induced events (on vertical axis) 
in six different arbitrary slices (red, green and blue columns) around the principal fault after 
different intervals of time (on horizontal axis). The distances of the limits of the slices are 
computed from the perimeter of the source and the reference magnitude in the simulation is 
6.0. Histogram (a) rapresents the distribution nearest to the source at shorter intervals of time 
(the time scale is purely qualitative). Histogram (b) rapresents the distribution farther from 
the fault at longer intervals of time. 



Figure (11): Spatial distribution of the density of induced events at t=0.01 year. 



Figure (12): Spatial distribution of the density of induced events at t=l year. 
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Figure (13): Spatial distribution of the density of induced events at t=100 year. 

Figure (14): Comparison between the temporal distribution of the cumulative number of 
induced events corresponding to synthetic data and the best fits obtained by the Omori law for 
the closer slices from the source. 

Figure (15): Comparison between the temporal distribution of the cumulative number of 
induced events corresponding to synthetic data and the best fits obtained by the Omori law for 
the farther slices from the source. 
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